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We report simulations of submonolayer epitaxial growth using a continuum phase field model. 
The island density and the island size distribution both show scaling behavior. When the capillary 
length is small, the island size distribution is consistent with irreversible aggregation kinetics. As the 
capillary length increases, the island size distribution reflects the effects of reversible aggregation. 
These results are in quantitative agreement with other simulation methods and with experiments. 
However, the scaling of the island total density does not agree with known results. The reasons are 
traced to the mechanisms of island nucleation and aggregation in the phase field model. 

PACS numbers: 68.35.Fx, 81.10.Aj, 81.15.Aa 



I. INTRODUCTION 

Epitaxial growth is an important phenomena that 
has attracted theoretical attention from many different 
points of view. The main motivation is to understand 
and predict the surface morphology as deposition pro- 
ceeds. Some calculations focus on the energy parameters 
that control individual adatom motion. [l|, |2fl Other cal- 
culations focus on the kinetic roughening of the surface 
that occurs after thousands of layers have been deposited. 
Q The sub-monolayer regime is particularly interesting 
because (i) comparison between experiment and theory 
can be used to extract diffusion and adatom detachment 
barriers and (ii) the kinetics of submonolayer growth is 
replicated in the subsequent multilayer regime. [1] 

Several theoretical methods have been used to study 
the kinetics of sub-monolayer epitaxial growth. The old- 
est of these exploit rate equations to predict total island 
densities and the distribution of island sizes in a mean 
field theory. [|, H Kinetic Monte Carlo (KMC) simula- 
tions are particularly popular because they are atomistic, 
they provide a visualization of the growing surface, and 
they make predictions that often agree with experiment. 
[2H9j A desire to avoid the computation-time restrictions 
of atomistic simulations led to the development of the 
continuum level set method (LSM), which focuses exclu- 
sively on the motion of steps. [13, [H| Level set simula- 
tions have been shown to reproduce the results of KMC 
simulations for both sub- monolayer total island densities 
and island size distributions. [H, [l3| 

A recent paper by Yu and Liu (TJ approached the sub- 
monolayer problem using a phase field method. Phase 
field modeling is a continuum approach to the kinetics of 
phase transformations which makes no use of atomistic 
information. For that reason, it is widely used to study 
evolution phenomena over large length and time scales 



that are inaccessible to other methods. \\Jj When ap- 
plied to the problem of step flow growth in the limit of a 
thin interface (between the solid and its vapor), the phase 
field model reduces to the classic step flow model of Bur- 
ton, Cabrera, and Frank. [l6| Yu and Liu wrote down a 
phase field model to study the density of islands in the 
sub-monolayer regime. They reported that this quantity 
scaled with the deposition flux F and the adatom surface 
diffusion constant D as N oc (F/D) 1 / 3 . This is the ex- 
pected result in the irreversible aggregation regime where 
island nucleate when two atoms collide and there is no 
detachment of atoms from island edges. 

The original motivation for this pa per was to repro- 
duce the island density results of Ref. [lj and to extend 
them to study the distribution of island sizes in the sub- 
monolayer regime. It turned out that our results differed 
from theirs in a interesting way which, we believe, demon- 
strates some of the virtues and some of the defects of the 
phase field method applied to this particular problem. 
Our main result is that the island size distribution shows 
scaling behavior. When the capillary length is small, the 
island size distribution is consistent with irreversible ag- 
gregation kinetics. As the capillary length increases, the 
islands size distribution reflects the effects of reversible 
aggregation. The results agree quantitatively with KMC 
and LSM simulations and with experimental data. The 
total island density scales with D/F, but the exponent 
is not |, nor does it change when the scaled island size 
distribution changes shape. 
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II. CALCULATION AL METHOD 



The phase field model of Yu and Liu uses two dimen- 
sionless variables, the adatom concentration u and the 
order parameter (surface profile) <p. These are coupled 
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by the evolution equations: 
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In (JlJ, the first term models the surface diffusion of 
adatoms. The second term models mass exchange be- 
tween the adatom population and the steps. The third 
term is the mean deposition rate, and the last term is 
a random variable which determines the points on the 
surface where deposited atoms land. In the term 
2 sin(27r0) identifies the terraces of the step profile with 
integer values of <fi. The term W 2 V 2 cf> determines the 
width W of the step which connects adjacent terraces and 
the term proportional to u — u eq causes the boundary of 
an island to move by the capture or release of adatoms. 
The final term in is a rate equation estimate of the 
island nucleation rate. 

To discuss our choice of parameters, we recall the 
"thin-interface" limit of the phase field model. [l7| This 
limit defines a capillary length and a kinetic coefficient /3 
from 
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where a\ — 0.36 and 02 = 0.51. More importantly, do 
and /3 are related to each other in exactly the same way 
as they are related in the Burton, Cabrera, and Frank 
model of step flow growth. 18] Namely, 



v = D[n • Vit] step = /3 1 [u-u eq - d n] 



step 



(5) 



where v is the velocity of a step, n is a unit vector nor- 
mal to the step, u eq is the equilibrium concentration of 
adatoms at a straight step, and n is the step curvature. 
The subscript "step" in ([5J means that the quantities in 
brackets are evaluated at the step edge. We consider the 
limit (3 — only, which corresponds to fast attachment of 
adatoms to step edges (surface diffusion limited growth). 
In that case, we get the Gibbs-Thomson equation Q 

[u]step = U eq + d [n} step , (6) 

and there is no loss of generality if we set u eq = 0. In the 
same /3 = limit, 



A = 



d 



and 



aia 2 M /3 
d Q D 
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In practice, we chose a unit length a and fixed W = a 
and D = 10 4 a 2 /sec. The free parameters of the model 



are do (units of a), F (units of ML/sec), and A„. We dis- 
cretized the coupled equations ((T|) and © on a L x L 
square lattice with L = 960 grid points and solved 
them using no-flux boundary condition at the lattice 
edges and a two-dimensional forward-time, central space 
(finite-difference) algorithm. A parallel algorithm (do- 
main decomposition) was used to speed up the compu- 
tation. We found good convergence using a spatial grid 
size Ax = 0.4a. The time step At is chosen so that 
At <C (Ax) 2 /D. To model depositions, we choose a grid 
site at random and set u = a 2 /(Ax) 2 at that site. We 
then repeated this step every 1/ '(FLAx) 2 seconds. The 
surface coverage is defined as 9 — Ft. 



III. RESULTS 
A. Nucleation & Aggregation 

Figure [1] illustrates the nucleation and aggregation be- 
havior produced by the phase field equations ([1]) and 
@. The left column shows the adatom density u at 
three successive times. The right column shows the or- 
der parameter <f> (surface morphology) at the same three 
times. Panel (a) shows the rapid, isotropic diffusion of 
the adatom concentration away from a deposition event 
which occurred at the point labelled (4). Through the 
nucleation term in ([2J , this distribution of u triggers the 
growth of a small spike in (j> at exactly the point (4). 
This spike, which we call a proto-island, is not yet visi- 
ble in panel (b) , which instead shows three proto-islands 
[labelled (l)-(3)] which were triggered by three earlier 
deposition events. The adatom density associated with 
these earlier events has completely diffused away by the 
time of deposition event (4). 

Understanding the fate of proto-islands is the key to 
understanding the behavior of the model overall. Some 
proto-islands grow into true islands by the capture of 
adatom density from other deposition events. Other 
proto-islands disappear because not enough adatom den- 
sity is captured before <f> itself "diffuses" away due to the 
interface width term W in . Our choice of W produces 
well-defined islands with sharp edges. Diffusion along the 
island edges is naturally included by the surface free en- 
ergy minimization that leads to ([2j . In detail, we label as 
a proto-island every set of one or more nearest-neighbor 
connected grid sites where 4> > 0.05. If the value of 4> 
at each connected site is called <pk, we form the quantity 
s = (Ax j a) 2 ~^2 k 4>k for each proto-island and monitor its 
value as time goes on. If s — > 0, we say that this proto- 
island has disappeared; if s > 1 we say this proto-island 
has become a true island composed of s atoms. 

Panel (c) in Figure Q] shows the expected adatom con- 
centration very soon after a deposition event at the point 
labelled (8). More interesting is panel (d), which shows 
seven true islands. Islands (l)-(3) evolved from the proto- 
islands (l)-(3) in panel (b). Islands (4)-(7) were produced 
by deposition events that occurred in the time between 
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FIG. 1: Time evolution of the order parameter and corre- 
sponding adatom concentration. 9 is the surface coverage. 
Note that the color bar is varied to optimize the contrast. 
For all panels, D/F = 10 7 , d = 1.44 x 10" 6 , A n = 8.4 x 10" 3 
and L = 80a. The surface coverage: panel (a) and (b) 
9 = 2.7 x 10~ 4 , panel (c) and (d): 9 = 2.2 x 1(T 2 , panel 
(e) and (f): 9 = 2.8 x 10" 2 



panels (a) and (c). A short time later, panel (e) shows 
that the adatom density associated with deposition event 
(8) has diffused entirely away. However, no island (8) has 
been created in panel (f) because proto-island (8) disap- 
peared. It did not grow to a true island because the ex- 
isting islands captured all the available adatom density. 
In other words, the island density in this neighborhood 
of the surface has saturated and further deposition only 
causes the existing islands to grow. Indeed, the very dark 
regions of panel (e) can be regarded as "denuded" zones 
around each island. 

The foregoing shows that the nucleation of an island 
in the phase field model occurs quite differently than it 
does in, say, an atomistic KMC simulation. There, de- 
posited atoms diffuse on the surface until they collide to 
form a stable island somewhere away from the deposi- 
tion point of either atom. We have said that the phrase 
"irreversible growth" is used if this collision produces a 
stable island. We speak of "reversible growth" if a just- 
nucleated island can dissociate back into adatoms. That 
being said, the aggregation behavior of the phase field 
model seems quite similar to that seen in KMC and LSM 



simulations. We will see in a moment that this similarity 
(dissimilarity) of the nucleation (aggregation) process to 
other simulation results has consequences for the behav- 
ior of the distribution of island sizes and for the total 
island density. 

For later use, we draws particular attention to the level 
set method to simulate sub-monolayer epitaxial growth. 
In LSM simulations, islands are nucleated at random po- 
sitions on the surface using a rate-equation-like weight- 
ing factor proportional to the square of the adatom den- 
sity. [13] The adatom density itself evolves as dictated 
by a uniform deposition flux at every point and a diffu- 
sion equation with specified boundary conditions at the 
moving edges of existing islands. The method is very 
computer-time intensive, but as mentioned earlier, the 
total island density and the distribution of island sizes 
agree very well with KMC simulations and with experi- 
ment. 



B. Island Size Distribution 

The island size distribution n s is the number of islands 
composed of s atoms. If s av is the average island size, it 
is well-known that a plot of the scaled quantity n s s^ v /9 
versus s/s av will collapse onto a single curve data col- 
lected for different values of D/F. One particular 
curve is characteristic of irreversible aggregation and the 
shape of this curve varies smoothly as the degree of re- 
versibility is increased by changing, say, the pair-bond 
energy in a KMC simulation. [8j 

Fig. [2Ja) shows island size distributions obtained from 
our phase field simulations model at very low coverage for 
D/F = 10 5 — 10 7 and various choices of the model pa- 
rameters do and A„. Each data point of the same symbol 
represents the average of at least 20 simulations. The 
scaling curve we find agrees very well with irreversible 
KMC and LSM simulations and with low temperature 
experimental data collected for Fe/Fe(001). Data 
collapse onto a single curve generally required us to re- 
duce the value of do as we increased the value of D/F. 
Doing this (or changing A„) produced very different total 
island densities, even though the scaled island size distri- 
butions were the same. For example, the data associated 
with the symbols ▲ and • in Fig. Ufa) have island densi- 
ties that differ by 25%. Similar behavior occurs in LSM 
simulations when the boundary conditions at the island 
edges are changed slightly. [H| Based on Figure , we 
conclude that the details of the island nucleation process 
are not critical to the shape of the island distribution 
when irreversible growth occurs. What matters is the 
subsequent process of monolayer capture by existing is- 
lands. 

Fig. [2^b) and Fig. [2{c) show the effect on the island 
size distribution of progressively increasing the capillary 
constant do- The A data in these two figures correspond 
to the same choices of D/F, A„, and 9 used in Fig.[^a). 
The change in shape we find for the scaled island size 
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FIG. 3: The island density at a coverage of 9 = 0.1 depends 
on both do and A n . Simulations are done on a lattice with 
1920 x 1920 grid points. D/F = 10 6 . 



the BCF problem of adatom diffusion on terraces, this 
simultaneously reduces the gradient of the adatom den- 
sity at the step edge in the leftmost equation in and 
thus retards the growth speed of an island. The capillary 
constant do measures the strength of the Gibbs- Thomson 
effect, Q which is the driving force for adatom detach- 
ment from step edges in phase field modeling. 



FIG. 2: The crossover scaling of island size distribution. 
Experimental data (large circles) are replotted from [lj| 
for different temperatures, and KMC data (open symbols) 
from @]. (a) ■: D/F = 10 5 ,d = 1.44 x 10" 4 ,6I = 
0.06, A„ = 0.03; • and A: D/F = 10 6 ,d = 1-44 x 10~ 5 
and 2.43 x 10 _5 ,A„ = 0.06 and 0.1,0 = 0.05 - 0.1. ♦: 
D/F = 10 7 , d = 1-44 x 10" 6 , A n = 8.4 x 10" 3 , 9 = 0.01. (b) 
■ : D/F = 10 5 ,d = 1.44 x 10 _4 ,6» = 0.1, A„ = 0.03; • and 
A: D/F = 10 6 ,d = 1-44 x 10~ 5 and 4.0 x 10 _5 ,A„ = 0.012 
and 0.1, 9 = 0.05 - 0.1. (c) D/F = 10 6 . A and ▼: 
d = 1.0 x 10 _4 ,A„ = 0.1, 9 = 0.05 and 0.1; ► and <: 
d = 3.2 x 10" 4 , A„ = 1,9 = 0.05 and 0.1. 



distribution as do increases agrees quantitatively with the 
change in shape seen in reversible KMC simulations when 
the pair-bond energy is decreased or (equivalently) when 
the critical island size is increased. [9j Our results also 
agree with reversible LSM simulations. [2fj| 

The step velocity in reversible LSM simulations is cal- 
culated from 



D[n ■ Vu} s tep ~ V d et, 



(8) 



where the second term takes account of the detachment 
of atoms from island boundaries. Typically, Vdet is taken 
to be proportional to the density of island edge atoms. 
This may be contrasted with our ([S]), which shows that 
increasing do has the effect of raising the adatom density 
at islands edges (which is zero in LSM simulations). For 



C. Total Island Density 

We have pointed out (in connection with Fig [1} that 
nucleation is treated rather differently in the phase field 
model than in KMC or LSM simulations. To emphasize 
this point, Fig.|3]shows the total island density as a func- 
tion of do and A n for D/F = 10 6 . The decrease in island 
density with increasing do is striking, but not hard to un- 
derstand. Larger do increases the relative magnitude of 
the first two terms of the right hand side of Eq. © , which 
preserves the equilibrium state (i.e. <f> = or <fi = 1). 
Consequently, proto-islands hardly grow in the beginning 
(when <f> is close to zero) and many of them diffuse away. 
The island density increases as A n increases also. This 
parameter is the coefficient of the nucleation term in ([2J . 
Given the same surrounding adatom concentration, as 
one adatom is deposited, a larger X n triggers a larger 
change of the order parameter, which is more likely to 
survive and become an island. 

The foregoing may be compared with a rate equation 
analysis or an LSM simulation, where the nucleation rate 
is determined by a global average of the adatom concen- 
tration over the whole domain. Specifically, 



dN/dt = Dcti(u 2 



(9) 



where a\ is the (constant) capture number. In the stan- 
dard rate theory of irreversible aggregations, © leads 
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FIG. 4: The island density scaling vs. D/F with different do 
and \ n . d = 1.44 x 10~ 6 - 3.25 x 10~ 4 , A„ = 0.0084 - 1, and 
9 = 0.1. 



to a well-known scaling law for the total island density: 
N ~ (D/F)~ x with x = 1/3- This is also seen in irre- 
versible LSM and KMC simulations. However, the mech- 
anism implied by (|9|) is not truly captured by (JlJ and 
([2|). Instead, our phase field model uses \ n Du 2 as a lo- 
cal estimate of the nucleation rate. We remind the reader 
that, unlike other simulation methods, most islands grow 
out of the initial adatom depositions in the phase field 
method. Be that as it may, upon fixing dp. and A„ and 
changing only D/F, we found that the total island den- 
sity shows distinct scaling behavior. This is shown in 
Fig. |U The curves of different color correspond to differ- 
ent values of do and A„ over a wide range. The average 
value for the scaling exponent is x ~ 0.65. It is worth 
remarking that the island size distributions from differ- 
ent data points on the same curve in Fig. 0] usually do 
not collapse very well. This suggests that the degree of 
reversibility is not the same. 

We do not fully understand the scaling seen in FiglU al- 
though we presume a simple analytic theory exists which 
can reproduce the observed exponent. On the other 
hand, we can gain some insight by looking into the time 
evolution of the island density in more detail. Fig. [5] 
is a typical curve of N(t) obtained from a phase field 
simulation with D/F = 10 7 . By changing the model pa- 
rameters as described in Fig. [3j we can match the island 
density produced by a KMC simulation with the same 
value of D/F. However, there is a clear discrepancy in 
the nucleation rate: the island density approaches the 
steady state much faster in our simulations than in the 
KMC simulations. In fact, all of our phase field simula- 
tions show similar behavior. Since the island size distri- 
bution is a characteristic of the aggregation regime, this 
could explain why we can obtain the scaling of island 
size distribution at a much lower coverage than expected 
from KMC simulations (see Fig. [5]). The fact that most 
islands tend to form at an earlier time is undoubtedly 
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FIG. 5: Nucleation shuts off faster in phase field simulations. 
For D/F — 10 7 and the same island density in the steady 
state, the time evolution of island density in phase field sim- 
ulations (black squares) reaches the steady state much faster 
than in KMC simulations (dashed line, replotted from [ITj]'). 
d = 1.44 x 10" 6 , A n = 8.4 x 10" 3 . 



caused by the initial adatom depositions (see Fig. [T]). It 
follows that the nucleation rate in this phase field model 
decreases faster than what we expect from Eq. © , which 
results in a stronger dependence on u and thus changes 
the scaling of the island density. 



IV. CONCLUSION 

In summary, we have shown that phase field mod- 
eling of sub-monolayer epitaxial growth reproduces the 
scaled island size distributions seen in experiment and ob- 
tained from other high-quality simulation methods. The 
crossover from irreversible aggregation to reversible ag- 
gregation is driven by the magnitude of a capillary con- 
stant which enters the Gibbs- Thomson equation. This 
shows that diffusion-limited aggregation phenomena are 
well-captured by the model. |21| On the other hand, the 
scaling of the island density itself disagrees with exper- 
iment and with other simulation methods. This implies 
that our model does not treat nucleation as accurately 
as one would like. One simple solution is to abandon the 
term \ n Du 2 in ^ and use the level set method strat- 
egy to nucleate new islands. We suspect this will produce 
the correct total island density without changing the high 
quality already obtained for the island size distributions. 
This might be important, moving forward, because the 
phase field method is less computationally intensive than 
the LSM and is much easier to implement at larger spatial 
scales and for more complicated epitaxial growth situa- 
tions. 
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